22 February 2021 | R. Rosa
We addressed, in the previous post Time average bounds via Sum of Squares, the problem of estimating the asymptotic limit of time averages of quantities related to the solutions of a differential equation.
Here, the aim is to consider an example, namely the Van der Pol oscillator, and use two numerical methods to obtain those bounds: via time evolution of the system and via a convex semidefinite programming using Sum of Squares (SoS), both discussed in the previous post.
This example is addressed in details in Fantuzzi, Goluskin, Huang, and Chernyshenko (2016). My main motivation is to visualize the auxiliary function that yields the optimal bound via SoS. That is the main reason to choose a two-dimensional system.
That bound depends on the chosen degree \(m\) for the auxiliary function appearing the SoS method. Here are the results for two values of \(m\):
We discuss, below, the details leading to these result and how to appreciate the plots above. (Notice you can rotate, pan and zoom the previous and subsequent images!)
The Van der Pol oscillator originated in the study of eletric circuits and appears in several other phenomena, from control, to biology and seismology.
In its simplest, and nondimensional, form, the equation reads
\[ {\displaystyle {d^{2}x \over dt^{2}} - \mu (1-x^{2}){dx \over dt}+x=0,} \]where \(\mu>0\). This equation has the stationary solution \(x(t)=0\), \(\forall t\in \mathbb{R}\), and all other solutions converge to a limit cycle that oscillates around the origin. The form of the limit cycle and the convergence rate to the limit cycle changes according to the value of the parameter \(\mu\).
Below you will find the dynamics of both \(x=x(t)\) and its derivative \(x'=dx(t)/dt\) for different values of \(\mu\), and with the initial condition \(x(0)=4.0\) and \(x'(0)=6.0\). Notice that, for small \(\mu\), the solution converges faster to the limit cycle, which is nearly a sinusoidal wave; while for large \(\mu\), the solution converges slightly slower and the solution develops spikes, as if firing up some signal information (think neurons in biological applications).
Another common representation of the dynamics of an autonomous system is that of phase portrait, in which we draw the orbit, or trajectory, in the phase space of the system. In this case, the phase space is \(xy\), where \(y=dx/dt\). We rewrite the second-order differential equation as a system of first order equations
\[ \begin{cases} x' = y, \\ y' = μ(1 - x^2)y - x. \end{cases} \]The same trajectories above, in the two extreme values for \(\mu\), are seen below in phase space.
The post Time average bounds via Sum of Squares addressed the problem of estimating the time-average
\[ \bar\phi(u_0) = \limsup_{T\rightarrow} \frac{1}{T} \int_0^T \phi(u(t))\;\textrm{d}t, \]for the solutions \(u=u(t), t\geq 0\), \(u(0)=u_0\), of a differential equation
\[ \frac{\textrm{d} u}{\textrm{d}t} = F(u), \]where \(\phi:X\rightarrow \mathbb{R}\) is continuous; \(F:X\rightarrow X\) is some locally Lipschitz function acting in some finite-dimensional space \(X=\mathbb{R}^n\), \(n\in\mathbb{N}\); and assuming the solutions generate a continuous semigroup \(\{S(t)\}_{t\geq 0}\), where \(S(t)u_0=u(t)\).
We exemplify this here with the Van der Pol system (2) and with the quantity
\[\phi(x,y) = x^2 + y^2. \]The most direct way to numerically estimate the bound is via a numerical evolution of the system, for a sufficiently long time, and taking the corresponding time average. Since the system has a globally attracting limit cycle (except for the unstable fixed point at the origin), we may simply consider a single trajectory for the estimate.
Notice we integrated for a very long time to have a proper convergence of the time average. Alternatively, knowning that a fixed fraction of the time integral over the time period converges to zero, we could start integrating at a later time and avoid the initial large values, due to the transient behavior, and the initial bumps, due to the periodic spikes in the solution and its derivative:
\[ \lim_{T\rightarrow \infty} \frac{1}{T}\int_0^T \phi(u(t))\;\textrm{d}t = \lim_{T\rightarrow \infty} \frac{1}{T}\left(\int_0^{T_*} \phi(u(t))\;\textrm{d}t + \int_{T_*}^T \phi(u(t))\;\textrm{d}t\right) \\ = \lim_{T\rightarrow \infty} \frac{1}{T}\int_{T_*}^T \phi(u(t))\;\textrm{d}t = \lim_{T\rightarrow \infty} \frac{T-T_*}{T}\frac{1}{T-T_*}\int_{T_*}^T \phi(u(t))\;\textrm{d}t \\ = \lim_{T\rightarrow \infty} \frac{1}{T-T_*}\int_{T_*}^T \phi(u(t))\;\textrm{d}t = \lim_{T\rightarrow \infty} \frac{1}{T}\int_0^T \phi(u(T_*+t))\;\textrm{d}t. \]In general, however, we may not have such a clear understanding of the dynamics so as to start from specific points or at specific times.
We now estimate the upper bound via optimization. As we have seen in Time average bounds via Sum of Squares, since \(\phi\) and \(F\) are polynomials, an upper bound is given by
\[ \sup \bar\phi \leq \inf\{C; \; (C,V)\in \mathbb{R}\times\mathrm{P}_m(X), \;C−ϕ−F⋅∇V = \texttt{SoS}\}, \]where \(\mathrm{P}_m(X)\) denotes the set of real polynomials on \(X\) with degree at most \(m\),and where, in this case, \(X=\mathbb{R}^2\). We expect the bound to become sharper as we increase the degree \(m\). Recall, as discussed in the previous post, the degree has to be even, otherwise it has no chance of being nonnegative.
Below is the result of the estimate, for some choices of \(m\):
| degree | bound |
|---|---|
| 4 | 274.527 |
| 6 | 5.817 |
| 8 | 5.567 |
| 10 | 5.022 |
Here is a corresponding plot, but skipping \(m=4\), for scaling reasons:
Now we get to the point of visualizing the auxiliary function obtained for each degree. Looking at (7) and its feasibility condition
\[ C−ϕ−F⋅∇V = \texttt{SoS} \geq 0, \]we see that \(C\) is smaller, the greater \(-F\cdot\nabla V\), i.e. the "closer" \(F\) points to \(-\nabla V\), or, in other words, the faster the orbit descends along \(V\), whenever possible.
The best situation is in a gradient flow, but that is not always the case. It may not even be possible to have \(F\) descend along \(V\) all the time; just think of a periodic orbit, with a nontrivial auxiliary function \(V\), such as in our case. Nevertheless, the longer the orbit descends along \(V\), the better. With that in mind, observe the figures below. You can rotate and zoom the figure as you wish, to better observe the behavior of the limit cycle with respect to \(V\). Notice how the condition just described improves as the degree of \(V\) is allowed to increase.